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Preface 


The  purpose  of  this  study  was  to  determine  the  feasibility  of  using 
Fast  Fourier  Transforms  (FFT)  to  solve  Boundary  Value  Problems  (BVP). 

Since  many  Boundary  Value  Problems  are  solved  using  some  type  of  Finite 
Difference  Method,  it  was  felt  that  a  comparison  between  these  two 
methods  might  provide  an  insight  into  the  usefulness  of  the  Fast  Fourier 
Transform  in  solving  BVP. 

The  one  dimensional  BVP  was  studied  to  help  provide  a  basis  for 
understanding  the  two  dimensional  BVP.  The  two  dimensional  BVP  was 
studied  using  only  a  simple  case  where  the  boundary  conditions  were  zero, 
but  could  be  easily  extended  to  non-homogeneous  boundary  conditions. 

The  understanding  and  insight  of  the  FFT  I  gained  these  last  twelve 
weeks  has  been  extraordinary.  My  thanks  to  Dr.  N.  Pagano  of  the  Air 
Force  Weapons  Lab  for  sponsoring  this  thesis.  I  would  like  to  acknowledge 
Dr.  Kaplan  for  his  never  ending  support,  even  when  the  prospects  of 
getting  the  one  and  two  dimensional  FFT  computer  codes  working  was  some¬ 
times  questionable.  I  would  also  like  to  thank  Cpt.  Ric  Routh  and 
Cpt.  Jim  Helton  who  provided  greater  understanding  and  appreciation  for 
the  FFTs .  And  finally,  I  must  thank  my  wife,  Jeanine,  and  my  children, 
Jennifer,  Greg,  and  Derek,  who  supported  and  sustained  me  through  this 

period.  — . . . . 
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Abstract 

The  purpose  of  this  study  was  to  determine  the  feasibility  of  using 
Fast  Fourier  Transforms  (FFT)  to  solve  Boundary  Value  Problems  (BVP)  and 
then  compare  the  results  to  those  of  the  Finite  Difference  Method  (FDM). 
Variations  of  Poisson's  one  and  two  dimensional  equations  were  used  as  a 
vehicle  to  develop  the  FFT  method.  For  the  one  dimensional  BVP,  both 
homogeneous  and  non-homogeneous  Dirichlet  boundary  conditions  were  con¬ 
sidered.  In  the  one  dimensional  BVP  the  inhomogeneous  function,  F(x), 
was  also  varied.  The  two  dimensional  BVP,  only  one  inhomogeneous  function, 
F(x,y),  and  homogeneous  boundary  conditions  were  used.  The  one  dimensional 
model  was  used  as  a  basis  for  developing  the  two  dimensional  model. 

The  analytical  solution  of  each  problem  was  compared  to  the  numerical 
solution  of  the  FDM  and  the  FFT  method  at  varying  mesh  sizes.  The  compu¬ 
tational  time  of  the  FDM  and  the  FFT  method  were  also  compared. 

The  results  indicate  that  the  FFT  is  extremely  efficient  in  the  two 
dimensional  BVP  because  of  the  computer  storage  space  required  and  the 
computational  time  needed  to  solve  the  FFTs.  The  accuracy  of  the  FFT 
compares  favorably  to  the  FDM  and,  as  the  mesh  size  decreases,  becomes 


more  accurate  than  the  FDM. 


FEASIBILITY  AND  COMPARISON  INVESTIGATION  OF  THE  USE 
OF  THE  FAST  FOURIER  TRANSFORM  AND  FINITE  DIFFERENCE 
METHOD  FOR  NUMERICAL  SOLUTION  OF  BOUNDARY  VALUE  PROBLEMS 

I .  Introduction 

Background 

Many  science  and  engineering  problems  require  the  solution  of  one  or 
more  Boundary  Value  Problems  (BVP).  Many  of  these  BVPs  cannot  be  solved 
analytically  because  of  irregular  boundary  conditions  or  complex  geome¬ 
tries.  The  most  common  method  of  solving  these  type  of  BVPs  is  by  means 
of  the  Finite  Difference  Method  (FDM).  In  this  method  the  derivatives  of 
the  partial  differential  equation  are  approximated  by  use  of  Taylor  series 
expansion,  which  reduces  it  to  a  set  of  algebraic  equations  that  can  be 
solved  by  simultaneous  equations.  The  solution  of  these  simultaneous 
equations  may  be  found  with  the  aid  of  a  computer  by  using  direct  matrix 
inversion  techniques  which  require  N3  operations,  where  N  is  the  number 
of  simultaneous  equations  (14:6).  Hence,  other  methods  which  require 
fewer  operations  have  been  developed  to  solve  sets  of  simultaneous  linear 


equations.  These  methods  include  both  direct  and  iterative  techniques 
and  in  each  case  decrease  the  computational  time  (3:111;  6:417). 

One  of  the  techniques  used  to  solve  BVPs  involve  using  Discrete 
Fourier  Transforms  (DFT).  The  computational  time  for  a  DFT  is  only  on 
the  order  N2,  but  to  solve  a  one  dimensional  BVP  it  becomes  N3.  When  the 
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Fast  Fourier  Transform  (FFT)  algorithm  is  applied  to  the  DFT  the  computa¬ 
tional  time  is  decreased  to  an  order  of  NLOG2N  operations,  which  equates 
to  4NLOG2N  +  4N  to  solve  the  same  one  dimensional  BVP  (1:8;  11:215). 

This  decrease  in  computational  time  greatly  enhances  the  use  of  FFTs 
over  direct  FDMs  such  as  Gauss  Elimination.  The  use  of  the  DFT  in  conjunc¬ 
tion  with  the  FFT  will  be  referred  to  as  the  FFT  method  throughout  the 
remainder  of  this  study,  since  the  DFT  takes  advantage  of  the  FFT  algorithm 
for  increasing  its  computational  speed. 

Purpose 

This  thesis  topic  stems  from  an  article  in  the  May  1984  Physics 
Today  (5),  which  mentioned  the  use  of  FFTs  to  solve  BVPs  in  two  and  three 
dimensions.  For  a  three  dimensional  BVP  with  a  grid  or  mesh  of  100x100x100, 
the  matrix  becomes  very  large,  namely  106x  106 .  This  size  matrix, 
according  to  the  article,  can  easily  be  solved  using  FFT  techniques 
(5:56).  A  closer  investigation  of  this  method  revealed  a  lack  of  published 
information  concerning  the  FFT  technique.  It  was  found  that  the  FFT 
method  was  faster  and  required  less  computer  space  than  conventional  FDMs 
for  two  and  three  dimensions  (5:547;  13:710).  In  the  information  that 
was  published  only  BVPs  with  homogeneous  boundary  conditions  were  solved. 
There  is  an  obvious  lack  of  information  on  the  feasibility  of  using  FFTs 
to  solve  BVPs  in  one,  two,  and  three  dimensions.  Additionally,  there  is 
no  known  published  information  concerning  the  use  of  FFTs  to  solve  BVPs 
with  boundary  conditions  other  than  zero. 

This  study  investigated  the  feasibility  of  using  FFTs  to  solve  BVPs 
and  compared  this  method  with  the  Gauss  Elimination  method  in  terms  of 
both  computational  speed  and  accuracy. 


Scope 

The  study  will  consider  only  the  problem  of  the  one  and  two  dimensional 
Poisson's  Equation,  with  Dirichlet  Boundary  Conditions.  Both  zero  and 
nonzero  boundary  conditions  will  be  analyzed.  In  each  dimension  the 
Poisson's  equation  will  be  solved  analytically,  and  numerically  using 
Gauss  Elimination  and  FFTs.  The  accuracy  of  the  Gauss  Elimination  method 
and  the  FFT  method  will  be  compared  to  the  analytical  solution.  The  compu¬ 
tational  time  of  the  Gauss  Elimination  and  the  FFT  method  will  then  be 
compared. 

Plan  of  Development 

The  initial  approach  was  to  develop  a  FFT  computer  program  to  solve 
a  one  dimensional  Poisson  Equation  with  boundary  conditions  equal  to 
zero.  A  Gauss  Elimination  ( GE )  program  was  then  developed  to  solve  the 
same  one  dimensional  Poisson  Equation  wit.i  boundary  conditions  equal  to 
zero.  These  two  programs  were  then  expanded  to  accept  a  two  dimensional 
Poisson  equation  with  boundary  conditions  equal  to  zero.  Modifications 
were  then  made  to  account  for  boundary  conditions  other  than  zero. 

The  accuracy  of  the  numerical  solutions  were  then  compared  to  the 
analytical  solutions,  and  computational  times  of  the  GE  method  were  com¬ 
pared  to  the  computational  times  of  FFT  method.  Finally,  limitations  on 
the  FFT  method  and  feasibility  of  possible  directions  of  further  research 
were  discussed. 
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II.  Theor 


In  order  to  understand  the  use  of  FFTs  to  solve  boundary  value 
problems,  it  is  necessary  to  understand  some  theory  about  the  Fourier 
Series  and  the  Discrete  Fourier  Transform. 

Fourier  Series 

Any  function  expanded  as  a  series  of  eigenfunctions  is  defined  as  a 

Fourier  Series,  where  the  interval  of  orthogonality  is  ( 0<  x  < 2ir  ) 
(17:65).  Because  the  series  is  periodic  any  2-jt  length  can  be  used.  If 
the  interval  ( 0<  x  <2tt  )  is  replaced  by  ( -L<  x  <L )  ,  then  the  Fourier 

Series  can  be  defined  as 


CD 

U  ( x )  =  a  /  2  +  V  [a  C0S(nrrx/L)  + 
o  u  „  n 

n  =  1 


The  coefficients  a  ,  a  and  b  are 
on  n 


b  sinfnirx/L)] 
n 


(2.1) 


a 

o 


a 

n 


b 

n 


l/L  J  U(x)dx 


-L 


b 

■/« 


l/L  J  U ( x ) cos ( nffx/L ) dx  n=l,2,3.  .  . 
-L 


u 

j  J'  U(  x  )sin( 


l/L  /  U( x  )sin( nirx/L  )dx  n=l,2,3.  .  . 

ft 

-L 


(2.2) 


(2.3) 


(2.4  ) 


This  same  series  can  be  moved  along  any  interval  ( -L<  x  <L)  or  from 
( 0<  x  <2L )  as  long  as  the  interval  remains  2ir  (  7:283  ).  Thus  the  Fourier 
coefficients  can  be  rewritten  as 


aQ  =  2/L J  U ( x ) dx 
o 


a 

n 


b 

n 


■h 


2/L  J  U(  x  )cos(  nirx/L  )dx  n=l,2,3. 
o 


u 

j  y  U(x)sii 


2/L  /  U(x)sin(niTx/L  )dx  n=l,2,3.  .  . 
o 


(2.5) 


(2.6) 


(2.7) 


Discrete  Fourier  Series 

The  Fourier  Series  can  also  be  expressed  as  a  discrete  finite  series. 
The  derivation  will  not  be  shown  here,  but  can  be  found  in  Numerical 
Analysis  books  by  Richard  W.  Hamming  (7),  or  by  Robert  Vichnevetsky  (16). 
There  are  three  orthogonality  relationships,  though,  that  aid  in  the 
derivation  of  the  Discrete  Fourier  Series  and  understanding  of  the  Discrete 
Fourier  Transform  (7:284). 


2N-1 

l  cos{  (2ir/L)k(Lp/2N)}cos{  (  2TT/L  )m(  Lp/2N  )  } 

p  =  0 


0  k^m 
N  k=m^O 

2N  k=m=0  (2.8) 


2N-  1 

l  cos{  (27T/L)k(Lp/2N)}sin{  (  2tt/L  )m(  Lp/2n  ) }  =  0  (2.9) 

p  =  0 


2N-1 

l  sin{  (2ir/L)k(Lp/2N)  }sin{  (  2tt/L  )m(  Lp/2N  ) } 

p=0 


0  k?!m 
N  k=m^0 

2N  k=m=0  (2.10) 


Based  on  these  orthogonality  relationships  the  Discrete  Fourier  Series  can 


be  defined  as 


iKi 


Is 

M 


U(x)  =  a  /2  +  £  {a  cos(k27rx/L)  +  b  sin  (  k2-jrx/L  )  } 


(2.11) 


ww 


where 
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a 


2N- 1 

a^  =  1/N  l  U ( x  )cos  (  2ukx/li ) 

k  =  0 


2N-  1 

=  1/n  £  U(x)sin(27rkx/L) 

"  k  =  0 


(2.12) 


(2.13) 


The  coefficients  a  and  b  are  also  called  the  Discrete  Fourier 
k  k 

Transforms  (16:50).  To  see  why  this  is  true  it  is  necessary  to  review 
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the  Fourier  Integral,  then  relate  the  Integral  to  the  Discrete  Fourier 
Transform. 


Discrete  Fourier  Transform 

For  the  purposes  of  this  study  the  Fourier  Integral  will  not  be 
derived,  but  only  stated.  Several  books  contain  the  derivation  of  this 
Integral,  which  include  E.C.  Titchmarch's  book  on  Fourier  Integrals  (15) 
The  Fourier  Integral  is  defined  as 


I 

) 

»» 

b 


u  ( w )  = 


( t  )e”iwt  dt 


(2.14-a) 


0(w] 


U(t)  =  1/2TT  /  U(w)e  dw 


(2.  14-b ) 
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where  w  *  2nf  .  The  function  U(w)  is  called  the  Fourier  Transform  of 
U(t),  and  U(t)  is  the  inverse  Fourier  Transform  of  U(w).  Equation  (2.14-a) 


m 


Figure  1.  Function  U(t)  on  the  Interval  ( 0<  t  <t  ) 

—  —  N 

and  (2.14-b),  also  known  as  the  Fourier  Integral  Theorem,  can  then  be 

discretized  to  obtain  the  Discrete  Fourier  Transform. 

Let  U ( t )  be  defined  as  some  function  between  the  interval  ( 0<  t  <t  ) 

—  —  n 

as  outlined  in  Figure  1.  One  can  then  determine  the  Discrete  Fourier 
Transform  of  U(t),  which  is  defined  as  U(w). 

If  one  allows  U(t)  to  be  discrete  for  every  Uft^),  where 
k-0,1,2.  .  .  N  ,  and  applies  the  Fourier  Integral  Theorem,  while  letting 

"FT"  be  defined  as  a  "Fourier  Transform,"  the  following  is  obtained. 

00 

FT{ U ( t ) }  =  U(w)  =  J U( t )e~1Wt  dt  (2.15) 

-  00 


In  equation  (2.15)  w  is  defined  as  the  frequency  w=2n/T  ,  and  T  is 
defined  as  the  length  from  0  to  t  .  Because  U(t)=0  for  t<0  then 


U(w) 


00 

yu(t)e-iwt  dt 

o 


(2.16) 


By  letting  U(t)=0  for  t>t  ,  then 


By  discretizing  the  integral  the  following  is  obtained 


N 

U(w)  =  J  0 


( t )e_iWt  dt 


(2. 18-a) 


U(w)  »  l  U(t  )« 


(2. 18-b) 


where  t  =N A  t=T  ,  since  N  is  defined  as  the  number  of  intervals  between 
N 

zero  and  t  and  At  is  defined  as  the  interval.  This  means  At=T/N  and 
N 

tk=kAt=kT/N  ,  then 


U(w)  =  £  U(t  )e_1W(kT/N)  T/N 


(2. 19-a) 


U(w)  =  T/N  £  U(t  ) 


-iw( kT/N ) 


(2. 19-b) 


Since  w=2ir/T  ,  then  w=mAw=m2-jT/T  ,  where  m=l,2,3.  .  .  .N  ,  which 
leads  to 


U(mAw)  =  T/N  l  U(t  ) 


-i(  2ttiti/T)  (kT/N) 


(2.20-a) 


U(mAw)  =  T/N  l  U(t  )e 


-i2irmk/N 


(2.20-b) 


where  the  Discrete  Fourier  Transform  is  defined  as 


U(w  )  =  FT{U(tk)} 


1/N 


Nfu(tk)e-l27Tkm/N 


(2.21-a) 


If  one  lets  "ft-1"  be  defined  as  the  Inverse  Discrete  Fourier  Transform, 
then 


U(tk)  =  FT-1  {u(w  )} 


N-l 

l 

m=0 


U  ( w  )  e 
m 


i2irkm/N 


(2.21-b) 


The  Discrete  Fourier  Transform  in  equation  (2.21-a)  is  called  the  complex 

form  of  the  DFT.  The  DFT  can  also  be  expressed  in  terms  of  sines  and 

i  i  x 

cosines  by  using  Euler's  identity,  e  =  cos(jx)  +  isin(jx)  ,  and  is 
defined  as 

N-l 

U(w  )  =  a  /2  +  V  [a  cos(27ikm/N)  +  b,  sin(  2TTkm/N )  ]  (2.22) 

mo  /,  1  k  k 

k  =  1 

It  is  obvious,  then,  from  equation  (2.11)  that  a  ,  a  and  b  are  the  same 

ok  k 

Fourier  coefficients.  Thus,  these  Fourier  coefficients  are  also  called 
the  Discrete  Fourier  Transforms  (16:50). 

Sampling  and  Aliasing 

There  are  two  terms  that  are  synonymous  with  DFTs  that  need  an 
explanation.  Since  the  DFT  is  not  continuous,  a  discrete  finite  number 
of  points  must  be  determined  at  which  the  DFT  will  be  calculated.  These 


equidistant  points  are  called  the  “sample"  over  which  the  DFT  is  evaluated. 


For  example,  U(t),  in  Figure  1,  is  being  “sampled"  over  the  interval 


w 


’2' 


The  other  term  that  needs  explanation  is  aliasing.  If  At,  in 
Figure  1,  is  too  large,  the  Fourier  coefficients  of  the  higher  frequencies 
will  fold  into  the  coefficients  of  the  lower  frequencies.  For  example, 
if  a  sample  of  eight  pcmts  is  taken  over  an  interval,  the  coefficient  of 
the  first  harmonic  will  be  equal  to  the  coefficient  of  the  seventh  harmonic; 
the  coefficient  of  the  second  harmonic  will  be  equal  to  the  coefficient  of 
the  sixth  harmonic;  and  the  coefficient  of  the  third  and  fifth  harmonic 
will  be  equal.  To  avoid  this  problem  of  aliasing  the  Nyquist  sampling 
rate  is  used.  This  formula  is  1/  At  =  2f  .  Simply  stated,  this  says  if 
At  <  l/2f  ,  then  aliasing  will  occur;  if  t  >  l/2f  then  aliasing  will 
not  occur  (1:85).  (The  f  is  defined  as  the  highest  frequency  component 
of  the  Fourier  transform. ) 


Fast  Fourier  Transform 

The  Fast  Fourier  Transform  (FFT)  algorithm  takes  advantage  of  the 
symmetry  of  the  trigonometric  functions  in  the  Discrete  Fourier  Transform 
(DFT).  The  regrouping  of  the  equations  in  calculating  the  DFT  reduces 
the  number  of  computational  operations.  The  Discrete  Fourier  Transform 
requires  on  the  order  of  N3  operations  to  solve  a  one  dimensional  BVP, 
where  N  is  defined  as  the  number  of  discrete  data  points.  The  Fast 
Fourier  Transform  algorithm  is  on  the  order  of  4NLOG2N  +  4N  operations  to 
solve  the  same  BVP.  The  FFT  requires  on  the  order  of  NLOGjN  operations, 
but  in  calculating  a  one  or  two  dimensional  BVP  the  number  of  computations 
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W-Y.i 
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must  include  calculations  of  the  FFT,  the  calculations  of  the  inverse 
FFT,  and  any  computations  conducted  while  the  function  is  transformed. 
Additionally,  any  modifications  to  the  FFT  algorithm  (i.e.,  deletions  of 
the  real  or  imaginary  components  of  a  complex  array)  add  additional 
calculations  (11:215).  The  theory  behind  the  FFT  algorithm  will  not  be 
discussed  in  this  study.  E.  Oran  Brigham's  book  on  FFTs  (1)  contains 
both  an  intuitive  and  theoretical  development  of  the  algorithm  and  is 
recommended  to  the  reader  who  wishes  to  gain  an  indepth  understanding  on 
how  and  why  the  FFT  algorithm  works. 


I. 


Poisson ' s  Equation  in  One- Dimension 


The  first  problem  examined  in  this  study  is  Poisson's  equation  in 
one  dimension.  The  general  form  of  the  equation  is 


d2U( x ) 
dx2 


-F(x) 


(3.1) 


where  U(x)  is  the  unknown  function  to  be  determined  and  F(x)  is  a  known 

function.  Both  homogeneous  and  non  homogeneous  boundary  conditions  are 
examined  along  with  different  values  of  F(x).  The  boundary  conditions  and 
the  function  F(x)  are  separated  into  four  distinct  cases. 


Case 

1 

U(0)=0 

U( L ) =0 

F(  x ) 

=  40 

Case 

2 

U(0)=0 

U( L ) =0 

F(  x ) 

=  X 

Case 

3 

U(  0 ) =2 

U( L ) =8 

F(  x ) 

=  40 

Case 

4 

U(  0 ) =2 

U{ L ) =8 

F(  x ) 

=  X 

Analytical  Solution 

The  general  solution  to  equation  (3.1)  for  F(x)  =  40  is  found  by 
direct  integration  and  takes  the  form 


U(x)  =  -20x2  +  Ax  +  B  (3.2) 

for  case  one  and  case  three.  Applying  the  boundary  conditions  in  case 
one,  equation  (3.2)  becomes 

U( x )  =  -20x2  +  200x  (3.3) 


and  applying  the  boundary  conditions  in  case  three,  equation  (3.2)  becomes 


U( x )  =  -20x2  +  200. 6x  ♦  2 


(3.4) 
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The  general  solution  to  equation  (3.1)  for  F(x) 


x  in  case  two  and  four 


is  also  found  by  direct  integration  and  takes  the  form 

U(x)  =  —  t  Cx  t  D  (3.5) 

D 


By  applying  the  boundary  conditions  in  case  two,  equation  (3.5)  becomes 


'  j  i  x  > 


-x3 

6 


1 0  0  x 
6 


(3.6) 


and  applying  the  boundary  conditions  in  case  four,  the  equation  becomes 


U  (  x ) 


259x 

15 


(3.7) 


Equations  (3.3),  (3.4),  (3.6),  and  (3.7)  are  the  analytical  solutions 
to  the  one  dimensional  Poisson  equation  (3.1). 


Numerical  Approximation 

This  same  equation  (3.1)  can  be  solved  using  numerical  approxima¬ 
tions.  Two  different  techniques  were  examined,  the  FDM  and  the  FFT 
method.  Both  methods  involve  the  subdividing  of  the  region  concerned  into 
N  nodes,  or  mesh  points,  and  solving  a  set  of  simultaneous  equations  for 
each  value  of  N.  The  FDM  can  solve  the  simultaneous  equations  using 
different  techniques  including  Gauss  Elimination,  Tridiagonal,  or  Iterative 
Methods.  The  FFT  method  uses  the  FFT  algorithm,  which  solves  the  trigo¬ 
nometric  equations  using  symmetry  of  the  sine  and  cosine  terms. 

In  general,  the  accuracy  of  the  solution  is  dependent  upon  the 
number  of  nodal  points  chosen:  the  finer  the  mesh,  the  larger  the  number 
of  nodes,  and  the  greater  the  accuracy. 


Figure  2.  General  Solution  to  Equation  (3.8)  with  4  Interior  Nodes 


I 

$ 


Finite  Difference  Method.  The  FDM,  as  mentioned  in  the  Introduction 
approximates  the  derivatives  of  the  differential  equation  by  use  of  a 
truncated  Taylor  Series  expansion.  The  derivation  can  be  found  in  many 
texts  which  include  C.F.  Gerald  and  P.0.  Wheatley  (6),  and  Clark  and 
Hansen  ( 3).  By  using  the  Taylor  Series  expansion  a  central  difference 
approximation  can  be  applied  to  equation  (3.1)  to  get 


d2U( x )  2U( x ) 

dx2  53  h  2 

X 


<Vi  '  2uk  *  Vi> 


-F(khx> 


(3.8) 


where  x  =  kh  .  The  h  is  defined  as  the  distance  between  nodal  points, 
x  x 

N  is  defined  as  the  total  number  of  interior  nodes,  and  k=l,2,3.  .  .N 

(14:6).  Equation  (3.8)  can  be  reduced  to  matrix  notation 


Au  =  -F 


(3-9) 


where  "A*  is  defined  as  the  coefficient  matrix  which  is  tridiagonal,  "u’ 
is  a  column  matrix  of  unknowns,  and  "F*  is  a  known  column  vector  (see 
Figure  2).  This  matrix  equation  can  then  be  solved  by  either  the  Tri¬ 
diagonal  method  (Thomas  method)  or  the  Gauss  Elimination  Method. 


M.  V 

v*v 


Fast  Fourier  Transform  Method.  The  FDM  takes  advantage  of  a  poly¬ 


nomial  interpolation,  whereas,  the  FFT  method  takes  advantage  of  a 


trigonometric  interpolation  (16:89).  In  the  FFT  method,  h=L/(N+l) 
where  L  is  defined  as  the  length  of  the  region  of  concern,  N  is  the 


number  of  sample  points  and  h  is  defined  as  the  length  of  each  sampled 


3 


interval.  It  follows,  that  x  =nh  ,  where  n=0,l,2,3.  .  .N+l  and  is 

n 

defined  as  the  set  of  grid  points  in  which  the  solution  of  equation  (3.1 


is  to  be  approximated.  By  inspection  the  BVP  must  be  an  odd  function. 


since  at  both  boundary  conditions  the  value  is  zero  or  a  constant.  The 
equation 


h 


U(x)  =  £  b  sin(k7TX/L)  k  =  1,2,3.  .  -N 


(3.10) 


can  be  used  to  approximate  the  solution  to  equation  (3.1).  It  follows 
from  equation  (2.11)  and  (2.13),  and  using  the  orthogonality  relationship 
of  equation  (2.10)  that 


b  =  2  /  (  N ♦  I  )  V  U  s  in  (  kirx/L  ) 
k  .  “ .  k 

<  =  1 


(3.11) 


By  taking  the  second  derivative  of  U(x),  equation  (3.10),  with  respect 


to  x  the  following  is  found. 


d2U( x ) 
dx2 


£  b  (  kn/L  ) 2  sin(  kiyx/L  ) 


(3.12) 


By  substituting  the  right  side  of  equation  (3.12)  into  equation  (3.1), 


m 


one  gets  the  following  results. 


N 

-  I  b(kTT/L)2sin(kTTx/L)  =  -F(x)  (3.13) 

k=l  k 


By  discretizing  x  to  x  ,  where  n=l,2,3.  .  .N  then  equation  (3.13) 

n 

becomes 

N 

V  b(k7T/L)2sin(kTrx./L)=F(x)  =  -f  (3.14) 

k  n  n  n 

k=  1 


Now,  one  can  solve  for  b^  using  the  orthogonality  relationship,  equation 
(2.10),  to  obtain 

N 

b  =  2/  (N+l)(L/krr)2  T  fsin(kTTx/L)  (3.15) 

k  L ,  n  n 

n  =  l 


where 


N 

FT(  f  )  =  f  =  2  /  ( N+ 1 )  7  f  sinfkTTx  /L)  (3.16) 

n  k  u .  n  n 

n=l 


Solve  for  U  by  discretizing  U(x)  such  that 
n 


U 

n 


N 


l 

k=  1 


b,  sin(k7Tx  /L ) 
k  n 


(3.17) 


The  specific  sequence  used  to  solve  for  each  nodal  value  of  using  the 
FFT  method  is  thus: 

1.  Compute  the  FT(f  ),  equation  (3.16),  by  using  the  FFT  algorithm. 

n 

Compute  each  value  of  b^  by  dividing  the  eigenvalue  of  U(x)  by 
the  FT ( f  ).  See  equation  (3.15). 


2. 


3.  Compute  the  nodal  value  of  U  by  using  the  inverse  FFT  algorithm, 

n 

equat ion  (3.17). 

Computer  Analysis 

FORTRAN  codes  were  developed  to  solve  the  one  dimensional  Poisson's 
equation  using  the  FDM  and  the  FFT  method.  The  codes  take  advantage  of 
subroutines  found  in  the  International  Mathematical  and  Statistical 
Libraries,  Inc.  library,  more  commonly  called  the  IMSL  library,  and  were 
run  on  the  Harris  800  main  frame  computer.  A  listing  of  the  codes  using 
the  FFT  method  are  in  Appendix  A.  The  FDM  used  both  the  Thomas  method 
and  the  Gauss  Elimination  method.  The  Thomas  method  was  programed  in 
FORTRAN  using  the  algorithm  found  in  Clark  and  Hensen,  page  47.  The 
Gauss  Elimination  method  was  programed  in  FORTRAN  using  the  IMSL  routine 
LEQIF. 

FFT  Algorithm.  The  IMSL  subroutine  used  throughout  this  study  of 
the  FFT  was  FFTCC .  The  FFT  algorithm  requires  a  complex  value  input  and 
provides  a  complex  value  output.  Additionally  it  computes  the  transforms 
over  a  2L  interval.  Because  of  these  peculiarities  it  was  necessary  to 
modify  the  computer  code  to  compute  the  FFT  for  the  one  dimensional 
Poisson's  equation.  Equation  (3.1)  is  described  only  over  an  interval 
of  L.  In  order  to  have  the  FFTCC  subroutine  operate  properly  over  this 
interval  it  was  necessary  to  input  the  function  over  a  2L  interval,  then 
only  accept  one  half  of  the  interval  in  the  output  as  the  solution.  For 
example,  in  case  one  the  function,  F(x)=40  ,  is  a  rectangular  function 

over  the  interval  0<  x  <_L  .  Since  equation  (3.1)  is  described  by  a  sine 


function,  which  is  odd,  it  was  necessary  to  input  F(x)=40  from  0<  x  <L 


and  F(x)=-40  from  K  x  <_2L  .  On  output,  only  the  REAL  value  of  the 


complex  ouput  was  accepted  from  0<  x  <L  .  In  the  FFTCC  subroutine  the 
first  coefficient  is  only  a^ .  In  electrical  engineering  terms  this  is 
referred  to  as  the  DC  component.  This  component  must  be  eliminated  when 
solving  the  one  dimensional  Poisson's  equation  since  the  solution  is  only 
a  series  of  sine  terms.  Hence,  in  the  main  FFT  program  it  was  necessary 
to  sum  for  b  from  n=2  to  N  rather  than  n  =  1  to  N. 

N 

b  =  y  [  2  /  ( N+ 1 )  ]  [  L  /  ( k-TT  >  1 2  f  sin(kTTX  /L)  (3.18) 

k  L  ,  n  n 

n  =  2 

Equation  (3.18)  eliminates  the  DC  component  and  allows  the  FFTCC  subroutine 
to  approximate  the  one  dimensional  Poisson's  equation. 

The  four  cases  outlined  on  page  12  vary  not  only  the  function  F(x), 
but  also  the  boundary  conditions.  Two  of  the  four  cases  used  boundary 
conditions  other  than  zero.  In  applying  non-zero  boundary  conditions  to 
the  FDM  the  non-zero  boundary  is  absorbed  into  the  F(x)  function.  This 
process  applies  similarly  to  the  FFT  method.  On  input  the  function  F(x) 
must  be  adjusted  at  the  endpoint  to  account  for  the  non-zero  boundary 
condition.  For  example,  in  case  three,  over  an  interval  from  0<  x  <L  , 
the  values  of  F(x)  would  be  as  listed  in  Table  I. 

By  programing  the  input  values  of  F(x)  in  this  manner  the  desired 
results  are  achieved.  This  process  seems  logical,  as  one  examines  a 
delta  function,  such  as  F(0)=2  ,  in  normal  space,  the  value  becomes  a 

constant  in  Fourier  space  and  bounds  the  equation  that  is  transformed. 

When  the  inverse  transform  is  conducted  this  constant  returns  to  a  delta 
function  adjusted  somewhat  by  the  computation  conducted  in  step  two  of 
the  FFT  method. 
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TABLE  I 


Values  of  F ( x )  vs  x  in  Case  3 
for  the  One-Dimensional  BVP 


Value  of  x 

Value  of  F(x) 

0 

2 

2 

40 

4 

40 

6 

40 

3 

40 

10 

8 

i 


i 


Exact  Solution.  Particular  solutions  to  equation  (3.1)  were  found 
by  letting  L=10  and  solving  for  x  at  values  of  2,  4,  6,  and  8.  The 
table  of  these  results  can  be  found  in  Appendix  B.  These  results  were 
then  compared  to  the  numerical  solution  using  the  FDM  and  the  FFT  method 
at  the  same  values  of  x,  2,  4,  6,  and  8. 

Nodal  Points .  Each  case  of  boundary  conditions,  as  described  on 
page  12,  was  computed  using  four  different  mesh  sizes.  The  mesh  size 
was  decreased  in  each  of  the  four  cases,  thus  resulting  in  a  better 
approximation  to  the  analytical  solution.  The  interior  nodal  mesh  points 
used  were  4,  9,  49,  and  99.  Appendix  B  contains  the  results  of  these 
computations . 

Average  Error .  The  accuracy  of  the  approximations  was  determined 
using  an  average  percent  error,  or  relative  error  (6:42).  For  the  one 
dimensional  case  the  average  percent  error  was  defined  as 


<E> 


Exact  Solution  -  Numerical  Solution 
Exact  Solution 


x  100  (3.19) 


average  percent  error  tor  each  set  of  interior  nodes  are  found  in 
Figures  3-6.  In  each  of  the  four  cases,  as  outlined  on  page  12,  the  FDM 
solution  was  found  to  be  the  same  as  the  analytical  solution.  Consequently 
the  only  error  differences  reflected  on  Figures  3-6  are  in  the  FFT  method. 
In  all  cases,  as  the  number  of  interior  nodes  increase  the  average  percent 
error  in  the  FFT  method  decreases.  The  computations  for  each  nodal  point 
for  the  FDMs  and  the  FFT  method  are  in  Appendix  B,  Tables  III-VI. 

Computational  Time.  In  order  to  calculate  N  interior  nodal  points 
using  the  FFT  method  it  is  necessary  to  calculate  the  Fourier  Transform 
of  2N  +  2  points.  This  is  due  to  the  interval  having  to  be  sampled  over 
the  entire  2L  period  as  described  in  the  paragraph  on  the  FFT  algorithm. 

A  comparison  of  computing  times  is  found  in  Figure  7.  It  is  important 
to  keep  in  mind  that  the  FFT  times  reflect  2N+2  calculations  while  the 
FDMs  are  only  N  calculations.  The  Thomas  method,  or  tridiagonal  method, 
proves  to  be  the  most  efficient  method  in  the  one  dimensional  problem. 

The  FFT  method  becomes  more  efficient  as  the  number  of  nodal  points  are 
increased  as  compared  to  the  Gauss  Elimination  ( GE )  method.  This  is  due 
to  the  4NLOGjN  +  4N  calculations  required  for  the  FFT  method  in  one  dimen¬ 
sion,  while  the  GE  method  requires  N3  calculations.  Figure  8  is  a 
comparison  of  the  total  computer  time,  which  includes  loading  of  arrays 
for  input  and  writing  the  solutions  to  a  file  on  output.  The  computation 
times  for  the  algorithm  and  the  computation  times  for  the  total  computer 
program  are  listed  in  Appendix  B,  Tables  VII  and  VIII  respectively. 


Figure  4.  Case  Two,  Average  Error  Comparison  Between 
FDM,  FFT  Method  and  the  Analytical  Solution 


'» 
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BOUNDARY  CONDITIONS  U(0)=2  &  U(10)=8 

£  FUNCTION  F(X)  •  X 


Figure  6.  Case  Four,  Average  Error  Comparison  Between 
FDM,  FFT  Method  and  the  Analytical  Solution 


Figure  8.  Comparison  of  Total  Computer  Time  for 

the  Thomas  Method,  the  GE  Method,  and  the 
FFT  Method 


Poisson* s  Equation  in  Two  Dimensions 


The  second  problem  examined  in  this  study  is  Poisson's  equation  in 
two  dimensions.  The  form  of  the  equation  is 


V*U(x,y)  =  -2 

with  boundary  conditions 
U( 0, y  )  =  0 

U(8.y)  =  0 

U(x, 0 )  =  0 

U(  x , 6  )  =  0 


(4.1) 

i  4 . 2-a ) 
(4.2-b) 
(4.2-c) 
(4.2-d) 


where  U(x,y)  is  the  unknown  function  to  be  determined.  This  BVP  is 
described  by  equations  (4.1)  and  (4.2)  as  a  rectangle  with  dimensions  of 
8  in  the  x-direction  and  6  in  the  y-direction  (see  Figure  9). 


Figure  9.  2-D  BVP  with  6  Interior  Nodes 
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TABLE  II 

Exact  Solution  to  2-D  Poisson's  Equation 
at  6  Interior  Nodes 


Nodal  Point 

Exact  Solution 

uu 

4.18250 

°12 

5.95685 

U13 

4 .  18250 

U21 

4.18250 

°22 

5.95685 

°23 

4.  18250 

Analytical  Solution 

Equation  (4.1)  can  be  solved  using  variable  substitution,  separation 
of  variables,  and  Fourier  Series  techniques.  The  complete  solution  to 
this  BVP  can  be  found  in  Appendix  C.  The  solution  to  equation  (4.1), 
after  applying  the  boundary  conditions  in  equation  (4.2)  is 

U (  x  ,  y  )  =  (  8x  -  x2  -  5  1 2 /tt  3  ) 


[sinh[(2n-l)TTy/8]  *sinh[  (2n-l)ir(6-y)/8] 
<2n-l)3sinh[3(2n-l)7T/4] 

sin{  (  2n- 1 ) ttx / 8  ] 


(4.3) 


The  exact  solutions  for  U(x,y)  at  the  6  interior  nodes  shown  in  Figure  9 
are  found  by  solving  equation  (4.3)  at  each  nodal  point.  A  simple  BASIC 
program  was  written  to  solve  for  these  nodal  points.  The  solutions  con¬ 
verged  to  the  fifth  decimal  place  after  only  six  terms  were  summed.  A 
summary  of  the  exact  solutions  are  in  Table  II. 
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Numerical  Solution 

The  equation  (4.1)  can  also  be  found  using  numerical  approximations 


as  outlined  in  Chapter  II.  The  mesh  is  superimposed  over  the  rectangle 
as  shown  in  Figure  9.  There  are  n  equally  spaced  nodes  in  the  x-direction 
with  a  mesh-size  of  h^,  and  m  equally  spaced  nodes  in  the  y-direction. 


with  a  mesh-size  of  h  .  The  total  number  of  interior  nodes  is  equal  to 

y 

m  x  n,  which  is  equal  to  N.  The  size  of  h  and  h  will  dictate  the  number 

x  y 

of  n  and  m  points,  hence,  the  number  of  simultaneous  equations  necessary 

to  be  solved.  The  accuracy  of  the  solution  in  general  is  dependent  on 

the  value  of  n  and  m  as  in  the  one  dimensional  problem. 

Finite  Difference  Method.  By  allowing  h  =h  and  representing  them 
-  x  y 

by  h,  the  solution  with  FDM,  using  central  difference  approximations,  can 
be  simplified.  Equation  (4.1)  can  be  approximated  by 


+  2  =  0 


(4.4) 


since  h  =h  =h  ,  then 

x  y 


l/h-^U  ,  -2U  +U  +U  - 2U  +U  )+2  =  0 

j + 1 , k  ;,k  j-l,k  l,l+k  j,k  3,1-k 


which  can  be  simplified  to 


1/hMU.  +U  +  U.  +U.  -4U.  ]+2  =  0 

3+1, k  3-1, k  j.l+k  3,1-k  ],k 


(4.6) 


Equation  (4.6)  is  a  matrix  equation  which  can  be  solved  directly  by  the 
Gauss  Elimination  method.  Figure  10  is  an  example  of  equation  (4.6) 
with  6  interior  nodes  (  3  nodes  in  the  x-direction  and  2  nodes  in  the 
y-direction ) . 
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Figure  10.  2-D  Solution  to  Equation  (4.6)  with  6  Interior  Nodes 


Fast  Fourier  Transform  Method 

The  solution  to  equation  (4.1)  using  the  FFT  method  is  quite  similar 

to  the  method  described  in  Chapter  III  for  the  one  dimensional  problem 

except  an  additional  dimension  is  added.  By  allowing  U(x,y)  and 

F(x,y)=2  to  be  extended  as  odd,  2L-periodic  functions  in  both  x  and  y, 

then  equation  (4.1)  can  be  solved  using  numerical  approximations  in  three 

steps  (16:151).  The  derivation  of  each  of  these  steps  can  be  found  in 

Appendix  D.  The  first  step  is  to  compute  the  coefficients,  F  ,  by  the 

]  * 

use  of  the  Fast  Fourier  transformation  of 


Jk 


4  ( h  h  )  M  N 

- — ^ —  Y  sin(  jumAx/x  )  T  sin ( kimAy/y  )f 

x  y  L.  max  L .  max  mn 


V  S  max 

max  max  m=l  n=l 


(4.7) 


where  x  is  the  maximum  value  of  x  in  the  x-direction  and  y  is  the 
max  max 

maximum  value  of  y  in  the  y-direction.  The  value,  f  ,  is  defined  as  the 

mn 

discrete  value  of  F(x,y),  where  h  =x  /(M+l),  the  mesh  size  in  the 

x  max 
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x-direction,  and  h  =y  /(N+l)  ,  the  mesh  size  in  the  y-direction. 

y  max 


M  is  defined  as  the  number  of  mesh  points  in  the  x-direction  and  N  is 

defined  as  the  number  of  mesh  points  in  the  y-direction.  The  second  step 


is  to  compute  the  coefficients,  C  ,  by  applying 

1  * 


jk 


■jk 


l  ( 37T / x  ) 2  +  (krr/y  ) 2 1 
max  max 


(4.8) 


The  final  step  is  to  compute  the  nodal  values  of  the  function  U(x,y)  by 
applying  the  inverse  FFT  from 


M-l 

N-  1 

U(  x, 

’  y )  =  l 

l 

j  =  i 

k=  1 

Computer 

Analysis 

V  y  C  s in{  jttx/x  )sin(k7ry/y 
=  1  k=l  ]k 


max 


(4.9) 


The  FORTRAN  code  for  the  FFT  two  dimensional  problem  is  listed  in 
Appendix  A,  and  uses  the  IMSL  Routine  FFTCC .  The  Gauss  Elimination  codes 
were  developed  using  the  IMSL  Routine  LEQT1F.  Because  of  the  straight 
forward  nature  of  the  LEQT1F  subroutine  the  FORTRAN  codes  for  the  GE 
method  are  not  listed  m  an  appendix. 

FFT  Algorithm.  There  are  three  important  points  that  must  be  con¬ 
sidered  when  applying  the  two  dimensional  FFT  algorithm.  The  first  two 
points  are  the  same  as  outlined  in  Chapter  III  under  the  FFT  Algorithm 
subheading.  First,  the  sampled  interval  must  be  extended  over  2L  to 
account  for  the  FFTCC  subroutine  that  makes  computations  over  a  2L  interval 
Second,  the  DC  component  must  again  be  removed,  as  in  the  one  dimensional 


problem,  prior  to  applying  step  two  (solving  for  C  coefficients).  This 

J  * 


31 


is  accomplished  by  starting  to  sum  at  j=2  and  k  =  2  rather  than  j=l  and 


k.  =  1  .  Finally,  care  must  be  taken  as  to  which  FFT  algorithm  is  imple¬ 
mented.  Ideally,  an  algorithm  that  computes  sine  FFTs  and  inverse  sine 
FFTs  would  be  desired.  The  author  used  the  IMSL  routine  FFTCC,  for  the 
two  dimensional  problem,  which  computes  the  FFT  of  a  complex  array  in  one 
dimension.  This  was  necessary  because  a  sine  FFT  that  computed  both  the 
since  FFT  and  the  inverse  sine  FFT  was  not  available.  The  FFTCC  subroutine 
uses  a  complex  kernel,  which  means  that  both  the  real  and  the  imaginary 
part  of  the  function  is  computed.  In  order  to  compute  the  two  dimensional 
FFT  of  a  sine  function  it  was  necessary  to  write  a  subroutine  to  compute 
the  two  dimensional  FFT  and  eliminate  the  real  portion  of  each  computation. 
Appendix  E  contains  a  mathematical  explanation  as  to  the  reasoning  behind 
this  computation.  With  the  cosine  terms  removed  the  transform  becomes  a 
double  sine  series  as  shown  by  equations  (4.7)  and  (4.8).  The  number  of 
computations  required  to  calculate  the  two  dimensional  FFT  then  becomes 
8NL0G_,N  +  ION  (11:217). 

Exact  Solution .  The  exact  solution  to  equation  (4.1)  was  found  at 
six  points  as  illustrated  in  Figure  9.  The  value  of  each  nodal  point  is 

listed  in  Table  II.  These  exact  solutions  were  then  compared  to  the  solu¬ 
tions  obtained  from  using  the  FDM  (Gauss  Elimination)  and  the  FFT  method. 

Nodal  Points .  Equation  (4.1)  was  solved  numerically  with  three 
increasingly  larger  number  of  nodal  points.  The  first  mesh  size  included 

35  interior  nodes,  which  resulted  in  the  h  and  h  being  equal  to  1.0. 

x  y 

The  next  mesh  size  included  165  interior  nodes,  which  resulted  in  the  h 

x 

and  h  being  equal  to  0.5.  The  final  mesh  size  included  713  interior  nodes. 
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which  resulted  in  the  h^  and  h^  being  equal  to  0.25.  Appendix  B  contains 
the  results  of  these  computations. 

Average  Error.  The  accuracy  of  the  approximations  was  determined 
using  the  same  average  percent  error  equation  as  discussed  in  Chapter  III. 

<E>  =  Exact  Solution  -  Numerical  Solution  x  100  (3.19) 

Exact  Solution 

Comparison  of  Approximations  to  Exact  Solutions .  The  plots  of  the 
average  percent  error  to  each  set  of  interior  nodes  are  found  in  Figure  11. 
As  the  number  of  interior  nodes  increases  the  accuracy  of  the  FFT  method 
increases.  Note,  that  the  accuracy  of  the  FDM  decreases  as  the  number  of 
nodes  increases.  This  decrease  in  accuracy  is  due  to  roundoff  error 
(6:38).  Table  IX,  Appendix  B,  contains  the  values  of  the  analytical 
solutions  as  compared  to  the  GE  and  FFT  methods. 

Computational  Time ■  The  comparison  of  the  computational  speed  between 
the  GE  method  and  the  FFT  method  can  be  found  in  Tables  X  and  XI,  Appen¬ 
dix  B.  Table  X  contains  the  values  of  the  algorithm  computational  time 
required  to  compute  either  the  GE  or  the  FFT.  Table  XI  contains  the 
values  of  the  total  computer  time  used  to  run  each  program.  Figures  12 
and  13  provide  a  graphic  illustration  of  the  computational  times. 

Figure  12  is  the  comparison  of  the  algorithm  computations,  whereas 
Figure  13  is  the  comparison  of  the  total  computer  time.  It  is  obvious 
that  the  N3  operations  required  for  GE  method  really  become  significant 
when  the  interior  nodes  are  increased  to  165,  as  compared  to  the  FFT, 
8NLOG2N  ♦  ION. 
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Comparison  of  Algorithm  Computational  Times 
for  the  GE  and  FFT  Methods  in  the  Two 
Dimensional  BVP 
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(LOGARITHM  SCALE) 


COMPUTATIONAL  TIME  ON 
THE  HARRIS  800  COMPUTER 


Not  only  does  the  computational  time  increase  for  the  GE  method, 
but  the  computer  storage  space  increases  dramatical ly .  For  example,  to 
compute  the  713  interior  nodes  for  the  GE  method,  it  requires  a  matrix 
713x713  plus  two  one  dimensional  matrices  of  713.  This  equates  to  over 
500,000  memory  storage  locations  (10).  On  the  other  hand,  to  compute  the 
same  713  interior  nodes  using  the  FFT  method,  it  requires  two  75x82 
matrices,  one  120x160  matrix,  one  160  matrix,  and  one  complex  120x160 
matrix.  This  equates  to  just  over  50,000  storage  locations.  (The  FFT 
matrices  are  computing  4x713  interior  nodes,  because  the  FFT  method 
requires  a  2L  interval  in  both  the  x-direction  and  the  y-direction .  ) 

This  large  storage  requirement,  in  addition  to  the  computational  time, 
detract  from  the  efficiency  of  the  GE  method  and  enhance  the  FFT  method. 
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V.  Conclusions  and  Recommendations 


Conclusions 

The  results  from  the  one  dimensional  problem  show  that  the  Thomas 
method  (Tridiagonal  Method)  is  the  most  efficient  of  the  methods  examined 
because  it  provides  the  exact  solution,  and  the  computational  time 
requires  only  8N  operations  (3:48).  The  FFT  method,  though  faster  than 
the  GE  method  as  the  value  of  N  increases,  cannot  provide  the  exact  solu¬ 
tion.  The  GE  method,  as  does  the  Thomas  method,  provides  the  exact 

solution  with  only  4  nodal  points.  This  means  small  memory  storage 
requirements,  and  small  computer  time  usage  in  the  one  dimensional  BVP. 
The  exact  solution  using  the  FDM  is  because  the  FDM  is  based  on  a  poly¬ 
nomial  approximation,  whereas  the  FFT  is  based  on  a  trigonometric 
approximation . 

The  power  of  the  FFT  method  is  not  realized  until  the  Poisson's 
equation  is  analyzed  in  two  dimensions.  The  Thomas  method  cannot  be 
utilized  in  the  two  dimensional  problem,  because  a  tridiagonal  cannot  be 
formed  in  two  dimensions  when  the  matrix  is  full  (5:56).  The  Thomas 

method  is  used  in  the  Hockney  method  for  two  dimensions,  where  the  two 

dimensional  problem  is  broken  down  into  a  one  dimensional  FFT  and  a  one 


dimensional  Tridiagonal  (8:95;  16:151).  Additionally  the  Thomas  method 
can  be  used  in  two  dimensions  if  the  matrix  is  sparsely  populated  (5:56) 
In  the  case  of  a  fully  populated  matrix  the  comparison  of  direct  numeri¬ 
cal  methods  is  narrowed  to  only  the  GE  method  and  the  FFT  method. 

The  average  percent  error  of  the  FFT  method  in  the  two  dimensional 
BVP  decreases  as  the  number  of  nodal  points  is  increased,  as  shown  in 
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Figure  11.  The  average  percent  error  of  the  GE  method  increases  with 
increased  nodal  points  due  to  roundoff  error  (6:38).  The  best  percent 
error  with  the  GE  method  is  0.27%  at  35  nodal  points,  whereas  the  best 
percent  error  with  the  FFT  method  is  0.92%  at  713  nodal  points.  One  would 
reason  then  that  the  Ge  method  is  more  efficient  since  it  provides  greater 
accuracy  with  less  nodal  points.  This  is  true,  but  when  a  larger  set  of 
nodal  points  are  required,  which  is  often  the  case  in  engineering,  the 
value  of  the  GE  method  is  greatly  reduced.  This  is  apparent  when  the  com¬ 
putational  times  and  the  computer  storage  requirements  of  the  GE  method 
are  compared  to  the  FFT  method  (see  Figures  12  and  13).  The  FFT  method 
is  3  times  faster  with  35  nodal  points  than  the  GE  method,  55  times  faster 
with  165  nodal  points,  and  1175  times  faster  with  713  nodal  points  (see 
Tables  X  and  XI,  Appendix  B).  Additionally,  the  computer  storage  require¬ 
ment  for  the  GE  method  is  10  times  greater  than  the  FFT  method  as  discussed 
in  Chapter  IV.  Hence,  not  only  is  the  FFT  method  faster  than  the  GE 
method,  but  it  requires  less  computer  memory.  The  FFT  method  seems  to 
provide  the  most  efficient  two  dimensional  method  of  solving  Poisson’s 
equation  numerically  with  full  matrices.  The  ADI  method  is  faster  than 
the  FFT  in  two  dimensions,  but  is  only  good  for  sparse  matrices  (5:56). 

If  the  two  dimensional  BVP  is  expanded  to  three  dimensions,  then  the 
efficiency  of  the  FFT  method  could  prove  to  be  even  greater. 

In  the  one  dimensional  problem  four  cases  were  considered  that  varied 
both  the  value  of  F(x)  and  the  boundary  conditions.  One  question  of  great 
concern  is  whether  or  not  the  FFT  method  can  accept  boundary  conditions 
other  than  zero  (13:697;  16:89).  This  study  found  that  the  FFT  method  can 
be  used  with  other  than  zero  boundary  conditions  in  the  BVP.  This  procedure 


was  discussed  in  Chapter  III  under  the  Fast  Fourier  Transform  Method. 

The  two  dimensional  Poisson  equation  was  only  studied  with  boundary  con¬ 
ditions  equal  to  zero.  It  is  logical  to  assume  that  non-zero  boundary 
conditions  can  be  applied  to  the  two  dimensional  BVP  and  solved  using  the 
same  methods  employed  in  the  zero  boundary  condition  problem.  The  diffi¬ 
culty  associated  with  this  would  be  to  describe  the  input  function  as 
was  discussed  in  Chapter  III.  The  same  procedure  would  be  used  as  was 
done  in  the  one  dimensional  problem,  except  the  function  would  be  expanded 
to  two  dimensions. 

Recommendations 

This  study  has  only  touched  the  surface  on  how  the  FFT  method  can  be 
used  to  solve  BVPs .  Only  boundary  conditions  of  zero  and  constants  were 
used  in  this  study  to  establish  a  basis  for  the  validity  of  this  method. 

The  FFT  method  needs  to  be  expanded  to  non-zero  boundary  conditions  in 
two  dimensions.  (J.  Rosser,  from  the  University  of  Wisconsin,  has  done 
some  work  in  the  area  (12:38).)  And,  then  expanded  to  three  dimensions 
as  mentioned  by  Fox  and  Otto  (5:56).  Additionally,  there  has  been  limited 
work  in  the  use  of  Neumann  boundary  conditions  (12:41;  13:707).  Additional 

research  needs  to  be  done  on  how  extensive  FFTs  can  be  used  in  conjunction 
with  Neumann  conditions.  In  this  study  and  in  all  the  references  found  on 
the  FFT  method,  the  only  equation  studied  was  the  Poisson  equation.  Hence, 
there  is  a  question  as  to  the  ability  of  this  FFT  method  to  provide  numer¬ 
ical  solutions  for  equations  other  than  elliptic  equations.  Research 
needs  to  be  done  to  determine  the  limitations  of  FFTs  on  solving  BVP 


Appendix  A:  Computer  Codes 


This  Appendix  contains  the  listing  of  the  FORTRAN 
programs  used  in  the  one  and  two  dimensional  cases  for 
computing  the  FFT.  All  FFT  FORTRAN  programs  were  run  on 
the  Harris  800  computer,  using  the  IMSL  routine  FFTCC. 

The  following  are  the  titles  of  the  programs  in 
Appendix  A  and  their  function. 


FUNCTION 

Compute  one  dimensional  Poisson 
equation  with  F(x)=40  and  Boundary 
Conditions  equal  to  zero  using  the  FFT 
algorithm. 

Compute  one  dimensional  Poisson 
equation  with  F(x)=40  and  Boundary 
Conditions  U(0)-2  and  U(10)=8  using  the 
FFT  algorithm. 


FFT2  Compute  one  dimensional  Poisson 

equation  with  F(x)=40  and  Boundary 
Conditions  equal  to  zero  using  the  FFT 
algorithm. 

NFFT2  Compute  one  dimensional  Poisson 

equation  with  F(x)=40  and  Boundary 
Conditions  U(0)=2  and  U(10)=8  using  the 
FFT  algorithm. 

Compute  the  two  dimensional  Poisson 
equation  with  F(x,y)*2  and  Boundary 
Conditions  equal  to  zero  using  the  FFT 
algorithm. 


TITLE 

FFT 

NFFT 


FFT2D1 


♦PROGRAM  NAME  FFT  F(X)  *  40  (REC  FUNCTION) 

♦AUTHOR  TODD  R  JONES 

♦DATE  14  OCT  1985 

********************************************************* 

♦THIS  PROGRAM  WILL  PROVIDE  A  NUMERICAL  SOLUTION  TO  A  ♦ 
♦ONE  DIMENSIONAL  BVP  USING  FFT.  THE  BOUNDARIES  MUST  BE  ♦ 
♦ZERO  AT  BOTH  ENDS.  THE  FUNCTION  F(X)  MUST  BE  KNOWN  AND* 
♦PLACED  IN  THE  PROGRAM.  THE  FFT  SUBROUTINE  IS  AN  IMSL  * 
♦SUBROUTINE  CALLED  FFTCC  AND  WILL  TRANSFORM  300  VALUES  * 
♦OF  N.  ♦ 

********************************************************* 

* 

*  DECLARATION  OF  VARIABLES 

* 

INTEGER  N , IWK( 1050) 

REAL  WK( 1050) ,L(300) 

COMPLEX  A ( 300) 

PI  =  4 . *ATAN ( 1 . ) 

* 

♦  INPUT  NUMBER  OF  POINTS  TO  BE  TRANSFORMED 

* 

PRINT*. 'ENTER  NUMBER  OF  DATA  POINTS  TO  BE  TRANSFORMED 
READ* ,  N 

* 

*  INPUT  FUNTION  DESCRIPTION  F(X) 

* 

A ( 1 )  =  (0,0) 

A  ( N/2 )  =•  (0,0) 

A(N)  =  (0,0) 

A ( 2 )  -  (40,0) 

A ( N/ 2+ 1 )  =  (-40,0) 

DO  10  I  =  3.N/2-1 
A( I)  =  A ( 2 ) 

10  CONTINUE 

DO  20  I  =  N/2+1.N-1 
A ( I )  =  A( N/2+1 ) 

20  CONTINUE 
* 

*  START  TIMER 

* 

CALL  BTIME 
DO  30  I  -  1,N 

A ( I )  »  CONJG( A( I )  ) 

30  CONTINUE 

* 

*  USE  FFT  SUBROUTINE  FFTCC 

* 

CALL  FFTCC(A,N,IWK,WK) 

DO  40  I  -  1  , N 

A( I )  -  CONJG( A (I))*2/N 
40  CONTINUE 

DO  50  I  -  2 ,  N 

A ( I )  -  ((10/((I-1)*PI))**2)*A(I) 

50  CONTINUE 
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CALL  FFTCC(A,N,IWK,WK) 
CALL  ETIME 

********  STOP  TIMER  ******** 

* 

*  PRINT  VALUES  OF  U(N) 

* 

OPEN(l ,FILE=’DAT1FT1 ' ) 
DO  60  I  -  l.N 

L( I )  -  REAL( A( I  )  ) 

60  CONTINUE 

DO  70  I  =  1  ,N/2 

WRITE  (1,100)  L ( I ) 

100  FORMAT(  '  * , F8 . 3  ) 

70  CONTINUE 
CLOSE  (1) 

END 


♦PROGRAM  NAME  NFFT  F(X)  «  40  (REC  FUNCTION) 

♦AUTHOR  TODD  JONES 

♦DATE  1  NOV  1985 

*********************************************************** 

♦THIS  PROGRAM  WILL  PROVIDE  A  NUMERICAL  SOLUTION  TO  A  * 

♦ONE  DIMENSIONAL  BVP  USING  FFT.  THE  BOUNDARIES  MUST  BE  * 
♦U( 0) =2 ,  AND  U( 10)=8 .  THE  FUNCTION  F(X)  MUST  BE  KNOWN  AND* 
♦PLACED  IN  THE  PROGRAM.  THE  FFT  SUBROUTINE  IS  AN  IMSL  * 

♦SUBROUTINE  CALLED  FFTCC  AND  WILL  TRANSFORM  300  VALUES  * 

♦OF  N.  * 

*********************************************************** 
* 

*  DECLARATION  OF  VARIABLES 

* 

INTEGER  N , IWK( 1050) 

REAL  WK( 1050) ,L(300) 

COMPLEX  A ( 300) 

PI  *  4 . ♦ ATAN ( 1 . ) 

* 

*  INPUT  NUMBER  OF  POINTS  TO  BE  TRANSFORMED 

* 

PRINT*, 'ENTER  NUMBER  OF  DATA  POINTS  TO  BE  TRANSFORMED' 
READ*  ,  N 

* 

*  INPUT  FUNCTION  DESCRIPTION  F(X) 

* 

A ( 1 )  *  (2.0,0) 

A ( 2 )  -  (40.0,0) 

A ( N / 2  )  =  (8.0,0) 

A( N/2+1 )  =  (-2.0,0) 

A ( N/2+2 )  -  (-40.0,0) 

A ( N )  =  (-8.0,0) 

DO  10  I  =  3.N/2-1 
A  (  I  )  =  A  (  2  ) 

10  CONTINUE 

DO  20  I  =  N/2+2, N-l 
A ( I )  =  A(N/2+2) 

20  CONTINUE 
* 

*  START  TIMER 

* 

CALL  BTIME 
DO  30  I  -  1  ,  N 

A ( I )  -  CONJG( A (  I ) ) 

30  CONTINUE 

* 

*  USE  FFT  SUBROUTINE  FFTCC 

* 

CALL  FFTCC(A,N, IWK.WK) 

DO  40  I  -  1,N 

A ( I )  -  CONJG( A(I))*2/N 
40  CONTINUE 

DO  50  I  -  2 ,  N 


50  CONTINUE 

CALL  FFTCC(AtN,IWK,WK) 
CALL  ETIME 

********  STOP  TIMER  ******** 
* 

*  PRINT  VALUES  OF  U(N) 

* 

0PEN(1,FILE='DAT1NFT1') 
DO  60  I  -  l.N 

L( I )  =  REAL( A( I ) ) 
PRINT*, I,'  \L(I) 

60  CONTINUE 

DO  70  I  =  l.N/2 
WRITE(l.lOO)  L ( I ) 

100  FORMAT (  '  ’ , F8 . 3 ) 

70  CONTINUE 
CLOSE  (1) 


♦PROGRAM  NAME  FFT2  F(X)  =  X  (RAMP  FUNCTION) 

♦AUTHOR  TODD  R  JONES 

♦DATE  14  OCT  1985 

********************************************************* 

♦THIS  PROGRAM  WILL  PROVIDE  A  NUMERICAL  SOLUTION  TO  A  ♦ 
♦ONE  DIMENSIONAL  BVP  USING  FFT .  THE  BOUNDARIES  MUST  BE  * 
♦ZERO  AT  BOTH  ENDS.  THE  FUNCTION  F(X)  MUST  BE  KNOWN  AND* 
♦PLACED  IN  THE  PROGRAM.  THE  FFT  SUBROUTINE  IS  AN  IMSL  * 
♦SUBROUTINE  CALLED  FFTCC  AND  WILL  TRANSFORM  300  VALUES  * 
♦OF  N.  * 

********************************************************* 
* 

♦  DECLARATION  OF  VARIABLES 

* 


INTEGER  N , IWK ( 1050) 
REAL  WK(  1050) ,L(300) 
COMPLEX  A ( 300 ) 

PI  -  4 . *ATAN( 1 . ) 


*  INPUT  NUMBER  OF  POINTS  TO  BE  TRANSFORMED 

* 


PRINT* , ' ENTER  NUMBER  OF  DATA  POINTS  TO  BE  TRANSFORMED' 
READ* , N 


*  INPUT  FUNTION  DESCRIPTION  F(X) 


A(0)  =  (0,0) 

A ( 1 )  =  -10 

A ( N/ 2 )  =  (0,0) 

A(N)  -  (0,0) 

DO  10  I  =  2.N/2-1 

A ( I )  =  A ( I - 1 )  +  20. /N 
10  CONTINUE 

DO  20  I  =  N/2+1.N-1 
A ( I )  =  A ( I - 1 )  +  20. /N 
20  CONTINUE 


*  START  TIMER 

* 


CALL  BTIME 
DO  30  I  =  1  ,  N 

A ( I )  =  CONJG( A( I ) ) 
30  CONTINUE 


*  USE  FFT  SUBROUTINE  FFTCC 


CALL  FFTCC(A,N, IWK.WK) 

DO  40  I  -  1,N 

A ( I )  -  C0NJG(A(I))*2/N 
40  CONTINUE 

DO  50  I  -  2 ,  N 

A ( I )  -  ( ( 10/ ( (I— 1)*PI) )**2)*A( I) 
50  CONTINUE 

CALL  FFTCC(A,N,IWK,WK) 
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CALL  ETIME 

********  STOP  TIMER  ********* 

* 

*  PRINT  VALUES  OF  U(N) 

* 

DO  100  I  -  1 , N 
L ( I )  -  REAL( A( I ) ) 
PRINT*, L(I) 

100  CONTINUE 
END 
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♦PROGRAM  NAME  NFFT2  F(X)  =  X  (RAMP  FUNCTION) 

♦AUTHOR  TODD  JONES 

♦DATE  6  NOV  1985 

********************************************************** 
♦THIS  PROGRAM  WILL  PROVIDE  A  NUMERICAL  SOLUTION  TO  A  ♦ 

♦ONE  DIMENSIONAL  BVP  USING  FFT.  THE  BOUNDARIES  MUST  BE  ♦ 
*U(0)»2,  AND  U( 10)=8 .  THE  FUNCTION  F(X)  MUST  BE  KNOWN  AND* 
♦PLACED  IN  THE  PROGRAM.  THE  FFT  SUBROUTINE  IS  AN  IMSL  ♦ 

♦SUBROUTINE  CALLED  FFTCC  AND  WILL  TRANSFORM  300  VALUES  * 

♦OF  N.  * 

********************************************************** 

* 

♦  DECLARATION  OF  VARIABLES 

* 

INTEGER  N , IWK( 1050) 

REAL  WK( 1050) , L( 300) 

COMPLEX  A( 300) 

PI  *  4 . *ATAN(  1 .  ) 

* 

*  INPUT  NUMBER  OF  POINTS  TO  BE  TRANSFORMED 

* 

PRINT* , ' ENTER  NUMBER  OF  DATA  POINTS  TO  BE  TRANSFORMED' 
READ* ,  N 

PRINT*, 'ENTER  VALUE  FOR  STEP  INCREMENT  "Z"' 

READ* ,  Z 

PRINT* , ' ENTER  VALUE  FOR  A(N/2+l)' 

READ*,A(N/2+2) 

* 

*  INPUT  FUNTION  DESCRIPTION  F(X) 

* 

A ( 1 )  =•  (2.0,0) 

A ( N/ 2 )  =.  (8.0,0) 

A ( N )  -  (-2.0,0) 

A ( N/2+1 )  =  (-8.0,0) 

DO  10  I  =  2 , N / 2  —  1 
A ( I )  -  A(I-l)  +  Z 
10  CONTINUE 

DO  20  I  -  N/2+2.N-1 
A  ( I )  =  A  (  I  —  1 )  +  Z 

20  CONTINUE 

* 

*  START  TIMER 

* 

CALL  BTIME 
DO  30  I  *  1 , N 

A ( I )  -  CON JG( A (  I ) ) 

30  CONTINUE 

* 

*  USE  FFT  SUBROUTINE  FFTCC 

* 

CALL  FFTCC(A,N,IWK,WI) 

DO  40  I  -  1 ,  N 

A ( I )  -  CON JG( A(I))*2/N 
40  CONTINUE 
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DO  50  I  =  2 ,  N 

A ( I )  =  ((10/((I-1)*PI))**2)*A(I) 
50  CONTINUE 

CALL  FFTCC( A , N , IWK , WK ) 

CALL  ETIME 

********  STOP  TIMER  ********* 

* 

*  PRINT  VALUES  OF  U(N) 

* 

DO  100  I  =  l.N 

L( I )  -  REAL( A( I ) ) 

PRINT*, I,'  ' , L( I ) 

100  CONTINUE 

OPEN  (1,FILE='DAT1NF2' ) 

DO  60  I  -  l.N/2 
WRITE( 1 , 200)  L ( I ) 

200  FORMATC  '.F8.4) 

60  CONTINUE 
CLOSE  (1) 

END 


Wf 


♦PROGRAM  NAME  FFT2D1  F(X,Y)=  2  (REC  FUNCTION) 

♦AUTHOR  TODD  JONES 

♦DATE  23  OCT  85 

************************************************************ 

♦THIS  PROGRAM  WILL  PROVIDE  A  NUMERICAL  SOLUTION  TO  A  TWO  (2) 
♦DIMENSIONAL  BVP  (POISSON  EQUATION)  USING  FFTs .  THE 
♦BOUNDARIES  MUST  BE  ZERO  ON  ALL  SIDES.  THE  FFT 
♦SUBROUTINE  IS  AN  IMSL  SUBROUTINE  CALLED  FFTCC. 

************************************************************ 

* 

♦  DECLARATION  OF  VARIABLES 

* 

INTEGER  IWK(75,82),N,IA,IJOB 
REAL  RWK(75,82) ,L(120, 160) 

COMPLEX  A( 120, 160) ,  CWK(160) 

COMMON  IWK , RWK , CWK , L , A 
PI  =  4.0  *  ATAN( 1 . ) 

* 

♦  INPUT  NUMBER  OF  POINTS  TO  BE  TRANSFORMED 

* 


PRINT*, ' 
READ* , N 1 

ENTER 

NO.  OF 

POINTS 

TO 

BE 

TRANSFORMED-X 

DIR’ 

PRINT*, ' 

ENTER 

NO.  OF 

POINTS 

TO 

BE 

TRANSFORMED-Y 

DIR' 

READ* , N2 
I A  1  =  N 1 
IA2  =  N2 


* 

*  ENTER  COMPLEX  ARRAY  F(X) 

* 

*****  QUADRANT  I  ***** 

DO  10  I  -  2.N1/2 
DO  10  J  =  2.N2/2 

A ( I ,  J )  =.  (2. 0,0.0) 
10  CONTINUE 
*****  QUADRANT  II  ***** 

DO  15  I  =  N 1/2+1, Nl-1 
DO  15  J  -  2 , N2/2 

A( I , J)  -  (-2. 0,0.0) 
15  CONTINUE 
*****  QUADRANT  III  ***** 

DO  20  I  =  2.N1/2 
DO  20  J  *  N2/2+1 , N2-1 
A ( I , J )  *  (-2. 0,0.0) 
20  CONTINUE 
*****  QUADRANT  IV  ***** 

DO  25  I  -  Nl/2+1 , N 1 - 1 
DO  25  J  -  N2/2+1 , N2-1 
A( I , J)  -  (2. 0,0.0) 

25  CONTINUE 
* 

*  START  TIMER 

* 


CALL  BTIME 


*  SET  I JOB 

* 

I JOB  =  -1 

* 

*  USE  IMSL  ROUTINE  FFT3D 

* 

CALL  TWODIM(A,IAl ,IA2,N1 , N2 , I JOB , IWK , RWK , CWK ) 

DO  30  I  =  2 , N 1 
DO  30  J  *  2 ,N2 

A ( 1 1 J )  -  (4*A(I,J))/((PI*(J-1)/X)**2+(PI*(I-1)/Y)**2) 
30  CONTINUE 
I JOB  =  +1 

CALL  T W 0 D I M ( A , IA1 , IA2.N1 , N2 , I JOB , IWK , RWK , CWK ) 

CALL  ETIME 

********  STOP  TIMER  ******** 

DO  35  I  =  1 , N1 
DO  35  J  -  1 , N2 

L( I , J)  =  REAL( A( I , J)  ) 

35  CONTINUE 
* 

*  PRINT  VALUES  OF  NODAL  POINTS  TO  A  FILE 

* 

OPEN( 1 ,FILE='DFT2D1 ' ) 

DO  40  I  =  1.N1/2 
DO  40  J  =  1.N2/2 

WRITE(l.lOO)  I , J , L( I , J) 

100  FORMATC  1  ,I4,I4,F8.3) 

40  CONTINUE 
CLOSE  (1) 

* 

*  PRINT  ARRAYS  TO  SCREEN 

* 

DO  45  I  -  1.N1/2 
DO  45  J  =  1.N2/2 

PRINT  200,1, J,L(I,J) 

200  FORMATC  ’  ,  14 , 14  ,  F8 . 3  ) 

45  CONTINUE 
END 

************************************************************ 

************************************************************ 

*******************  SUBROUTINE  TWODIM  ********************** 
************************************************************ 

************************************************************ 

SUBROUTINE  TWODIM( A , I A 1 , 1 A2 , N 1 ,  N2 , 1 JOB , IWK , RWK , CWK  ) 
INTEGER  IA1,IA2,N1 , N2 , I JOB , I WK ( 1 ) 

REAL  RWK ( 1 ) 

COMPLEX  A ( I A 1 ,IA2) ,CWK( 1) 

INTEGER  I,J,K,L,M,N 
REAL  R 1 2 
COMPLEX  C 1 2 

* 

*  DETERMINE  TRANSFORM  OR  INVERSE  TRANSFORM 


*  INVERSE  TRANSFORM 

* 

DO  15  I  *  1,N1 
DO  15  J  «  1 , N2 

A( I , J )  -  CONJG(A(I, J)) 

15  CONTINUE 

* 

*  TRANSFORM  SECOND  SUBSCRIPT 

* 

10  DO  20  L  -  1 ,  N 1 

DO  25  M  =  1,N2 
CWK(M)  =  A ( L , M ) 

25  CONTINUE 

CALL  FFTCC( CWK , N2 , IWK , RWK ) 
DO  30  J  *  1 , N2 

A(L, J)  -  AIMAG(CWK( J) ) 

30  CONTINUE 

20  CONTINUE 
* 

*  TRANSFORM  FIRST  SUBSCRIPT 

* 

DO  35  J  -  1.N2/2 
DO  40  K  -  1 , N1 
CWK(K)  -  A ( K , J) 

40  CONTINUE 

CALL  FFTCC(CWK,N1,IWK, RWK ) 
DO  45  L  »  1 , N 1 

A(L , J)  *  AIMAG(CWK(L) ) 

45  CONTINUE 

35  CONTINUE 
* 

*  INVERSE  TRANSFORM 

* 

IF  ( IJOB.GT.O)  GO  TO  55 

R 1 2  -  N 1 *N2 
C 1 2  -  CMPLX( R 1 2 , 0 . 0) 

DO  50  I  -  1 , N 1 
DO  50  J  -  1,N2 

A ( I • J )  -  C0NJG(A(I,J))/C12 
50  CONTINUE 
55  RETURN 
END 


Appendix  B:  Tables  of  Data 

This  Appendix  contains  the  tables  of  data  collected 
from  the  one  dimensional  and  the  two  dimesional  Poisson's 
equation  to  include  Average  Errors  and  Computational  Times. 


Table  III 

Average  Error  for  the  One  Dimensional  Poisson  Equation 
with  F( x ) =40  and  U(0)=U(10)=0 


Analytical 


4-Interior  Nodes 


371.97 

468.98 
445.84 
302.58 


Total  Average  Error 


9-Interior  Nodes 


Analytical 


360.71 
483.90 
468.87 
320  I  315.64 


Total  Average  Error 


49-Interior  Nodes 


Average  Error  (Z 


7.77 


FFT 


72 

81 

32 


1.36 


4.30 


Analytical 


Average  Error 


*>1 


320 

480 

480 

320 

320 

480 

480 

320 

331.32 

482.85 

479.08 

320.02 

0 

0 

0 

o 

[Total  Average  Error! 

0  I 

1  1.08  i 

99-Interior  Nodes 


Analytical 


316.35 

479.19 

481.23 

322.47 


Average  Error 


FDM 


0 

0 

0 

0 


Total  Average  Error!  0 


0.58 


V/t 


m 


Table  IV 


loHiBmnm 

Ruiu:u 

FDM 

32 

32 

34.97 

0 

56 

56 

60.66 

0 

64 

64 

68.55 

0 

48 

48 

50.65 

0 

Total  Average  Error 

0 

mil 


F 


9.28 

8.32 

7.11 

5.52 


7.56 


49-Interior  Nodes 


X 

■mnnnmi 

2 

4 

32 

56 

32 

56 

6 

64 

64 

64.98 

8 

48 

48 

48.58 

n.m  w.uj  Jtumrni 


Average  Error  (%. 


FDM  I  FFT 


0 
0 
0 

0  I  1.21 


0  1.59 


32 

32 

32.3 

0 

56 

56 

56.49 

0 

64 

64 

64.49 

0 

48 

48 

48.29 

0 

Total  Average  Error 


0.78 


WlCM  sr  vO  00  I  XlCN  vO  COl  I  X  <N  <t  00  \ 


Table  V 


Average  Error  for  One  Dimensional  Poisson's  Equation 
with  F(x)=40  and  U(0)=2  and  U(10)=8 


4-Interior  Nodes 


Analytical 


323.2 

484.8 
485.6 

326.8 


323.2 

484.8 
485.6 

326.8 


349.56 

441,31 

373.51 

152.10 


e  E 


Average  Error  (Z 


FDM  FFT 


0 

0 

0 

0 


17.57 


9-Interior  Nodes 


Analytical!  F 

U 


323.2  323.2  352.78 

484.8  484.8  478.75 

485.6  485.6  444.76 

326.8  326.8  250.86 


Total  Average  Error 


FFT 


15 

25 

41 


23.23 


10.51 


49-Interior  Nodes 


Analytical 


323.2 

484.8 
485.6 

326.8 


323.2 

484.8 
485.6 

326.8 

305.35 

474.55 

483.75 

332.95 

0 

0 

0 

0 

Total  Average  Error! 

0 

Average  Error  (Z 


FFT 


5.52 

2.11 


0.38 

1.88 


2.47 


99-Interior  Nodes 


Analytical 


485.6 

326.8 


323.2 

484.8 
485.6 

326.8 


312.89 

477.49 

482.09 

326.69 


0.03 


1.36 
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Table  VI 

Average  Error  for  One  Dimensional  Poisson's  Equation 
with  F(x)»x  and  U(0)=2  and  U(10)=8 


4-Interior  Nodes 


nm 

■moi 

S33HI1 

FDM 

FFT 

35.2 

35.2 

42.46 

0 

20.62 

60.4 

60.4 

62.98 

0 

4.27 

69.6 

69.6 

70.56 

0 

1.38 

54.8 

54.8 

59.64 

0 

8.83 

Total  Average  Error 

0 

8.77 

9-Interior  Nodes 


mm 

SZERill 

nm 

FDM 

FFT  I 

35.2 

35.2 

34.56 

0 

mmm 

60.4 

60.4 

63.78 

0 

69.6 

69.6 

66.07 

0 

:  -  ‘  MM 

54.8 

54.8 

57.52 

0 

4.46 

Total  Average  Error 

0 

4.23 

49-Interior  No d e s 


mawm 

HQQSHHI 

■SRSH 

IBfSIV  |  ■ 

Average 

irror  (%) 

nn 

FDM 

FFT 

35.2 

0 

mmsm 

0 

69.6 

69.6 

68.05 

0 

54.8 

54.8 

54.53 

0 

Total  Average  Error 

0 

1.18 

99-Inter i or  N odes 


■m 

■KjQZHfil 

■Him 

■QjH 

Average  i 

irror  (I) 

■OOH 

FDM 

FFT 

35.2 

35.08 

0 

■ 

60.4 

60.39 

0 

1 

69.6 

69.6 

68.10 

0 

' 

54.8 

54.8 

54.03 

0 

■nos 

Total  Average  Error 

0 

0.97 

Table  IX 


Average  Error  for  the  Two  Dimensional  Poisson  Equation 
with  F(x)=2  and  all  Boundary  Conditions  equal  to  Zero 

35-INTERIOR  NODES  _ 


Nodal 

Analytical 

FDM 

FFT 

Ave  Error  (% 

Points 

U(x , y ) 

U(x,y) 

U(x, y ) 

FDM 

FFT 

Ull 

4.8125 

4.794 

4.742 

0.38 

1.46 

Ul  2 

4.8125 

4.794 

4.756 

0.38 

1  .17 

U21 

5.9568 

5.960 

6.111 

0.05 

2  .  58 

U22 

5.9568 

5.960 

6.126 

0.05 

2.84 

U31 

4.8125 

4.794 

4.625 

0.38 

3.89 

U32 

4.8125 

4.794 

4.972 

0.38 

3.31 

iTotal  Average  Error 

0.27 

2.54 

165-INTERIOR  NODES 


Nodal 

Analytical 

FDM 

FFT 

Ave  Error  (* 

Points 

U(x.y) 

U(x,y) 

u(*,y) 

FDM 

FFT 

Ull 

4.8125 

4.856 

4.705 

0.90 

9.19 

Ul  2 

4.8125 

4.856 

4.808 

0.90 

0.09 

U21 

5.9568 

6.026 

5.999 

1.16 

0.71 

U22 

5.9568 

6.026 

6.136 

1.16 

3.00 

U31 

4.8125 

4.856 

4.731 

0.90 

1.69 

U32 

4.8125 

4.856 

4.836 

0.90 

0.48 

ITotal  Average  Error 

0.98 

2.36 

713-INTERIOR  NODES 


Nodal 

Points 

Analytical 

U(x,y) 

FDM 

U(x,y) 

FFT 

U(x.y) 

Ave  Error  (% 
FDM  |  FFT 

Ull 

4.8125 

4.872 

4.771 

1.23 

0.86 

Ul  2 

4.8125 

4.872 

4.801 

1.23 

0.23 

U  2 1 

5.9568 

6.043 

6.036 

1.44 

1 . 32 

U22 

5.9568 

6.043 

6.074 

1.44 

1.96 

U31 

4.8125 

4.872 

4.770 

1.23 

0.88 

U32 

4.8125 

4.872 

4.799 

1.23 

0.28 

ITotal  Average  Error 

1.30 

0.92 
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Appendix  C:  Analytical  Solution  to  the  2-D  BVP 


The  solution  to  the  BVP  described  by  equation  (4.1)  and  (4.2)  is 
solved  by  using  a  variable  substitution  and  the  Fourier  Series  method  (4), 
By  making  a  variable  substitution 


U (  x ,  y )  =  V( x , y )  -  x 2  +  8 x 


( C . 1-a ) 


32u 

3x2 


32v 

3x2 


( C . 1-b ) 


32U 

3y2 


32v 
3y 2 


(C. 1-c) 


V(O.y) 


0 


(C. 5-a  ) 


V(  8, y )  =  0 


(C. 5-b) 


V( x, 0 )  =  x2  -  8x 


(C.5-C) 


V( x , 6  )  =  x2  -  8x 


(C. 5-d) 


By  using  separation  of  variables,  equation  (C.3)  leads  to  a  Sturm-Liouville 

problem  that  can  be  solved  using  the  Fourier  series.  The  solution  takes 
the  form 


V( x , y )  =  l 
n=  1 


A  sinh(niry/8)+B  sinh[  (nn/8 )  ( 6-y )  ]  .  , 

n _ n _  sin(nnx/8) 


sinh(  3nir/4  ) 


( C  .  6  ) 


This  can  be  verified  by  referring  to  Churchill  and  Brown,  page  136  (2:136) 
The  coefficients  A^  and  B^  can  be  solved  by  applying  boundary  conditions 
from  equation  (C.5). 


V(x,y)  =  Y  B  sin(nirx/8)  =  x2  +  8x 

L  i  n 

n=  1 


( C . 7-a ) 


V(x,y)  =  y  A  sin(n7ix/8)  =  x2  +  8x 

L .  n 
n=l 


( C. 7-b ) 


Thus,  A  =  B  and  can  be  solved  as  follows, 

n  n 


A^  =  2/8  J  (  x2 -8x  )sin(  nTjx/8  )dx 
0 


( C . 8-a ) 


An  -  1/4  /  x2sin(  njrx/8  )dx  -2  J  xsin(  nirx/8  )dx 


( C . 8-b ) 


A^  =  1/4  -  (  8x2 /ivjt  )cos  (  n7ix/8  )  +  12 8x/ (  n2TT2  ) sin (  mrx/8  ) 


+  1024 /  ( n3  it3  )cos  (  nyrx/8  ) 


JO 


-  (  8x/nir  )cos(  nirx/8 ) 


+  64/(n27T2  )sin( njrx/8 ) 


( C . 8-c  ) 


A  =  l/4[-512(-l)n/(nrr)  +  1024(-l)n/(n3Tt3  )  -  1024/ (  n3TT3  )  ] 
n 


-  2  [  -  (  64/niT )  (  - 1 )  ]  ( C .  8-d  ) 


An  =  Bn  =  256/(n3?73  )[  <  - 1 )"- 1  ]  (C.8-e) 

By  substituting  A^  and  into  equation  (C.6)  the  following  is  obtained. 

_  2561  ( - 1 ) n- 1 1  ?  t  sinh(  rnry/8 )+sinh{(mr/8)(6-y)}l 

,y  n3u3  sinh(3nn/4) 

n  =  1 

sin(nny/8)  (C.9) 

Equation  (4.1)  is  solved  by  making  the  variable  substitution  from  equation 
(C.9)  into  equation  (C.l-a). 

U(  x,  y)  =  (  8x  -  x2  -  5 1 2 / it 3  ) 


{sinh(  (  2n-l  )7ry/8  ]  »sinh[  (2n-l)7T(6-y)/8)} 
(2n-l)3sinh[3(2n-l)n/4] 


sin [  (  2n- 1  )ttx/8  1  (C.10) 


Equation  (C.10)  was  solved  by  programing  the  equation  into  BASIC  and 
running  the  program  on  a  Z150  PC.  It  was  found  that  the  solution  at 
specific  nodal  points  converged  after  six  iterations  to  the  fifth  decimal 
place.  The  specific  nodal  point  solutions  can  be  found  in  Table  II. 


Appendix  D:  Solution  to  the  2-D  FFT  BVP 


The  solution  to  equation  (4.1)  using  the  FFT  method  is  simply  an 
extension  of  the  one  dimensional  problem  described  in  Chapter  III.  By 
allowing  U(x,y)  and  F(x,y)  to  be  extended  as  odd,  2-periodic  functions  in 
both  x  and  y,  the  general  form  of  equation  (4.1)  is 

V  2  U  (  x ,  y  )  =  -  F  {  x .  y  :■  (  D .  1  ) 

and  can  be  represented  by  the  Fourier  series 

N-l 

I  F  s  in  (  jiix/x  )sin(kny/y  )  (D.2) 

lk  max  max 

N-l 

V  C  sin(j7ix/x  )  sin  (  kit  y/y  )  (D.3) 

lk  J  max  max 


where  C  ,  is  defined  as  the  coefficient  of  the  Fourier  series  and  the 
Jk 

mesh  size  h  =x  /(M+l)  and  h  =y  /(N+l)  (17:149).  The  boundary 

x  max  y  'max 

conditions  in  equation  (D.l)  are  zero  on  all  boundaries,  as  are  the 
boundary  conditions  of  equation  (4.1).  Now,  by  taking  the  second  deriva¬ 
tives  of  U(x,y)  with  respect  to  x  and  then  y  the  following  is  obtained. 


M-l 

F(x,y)  =  l 
1  =  1 


and 


M-l 

U(x.y)  =  l 
1  =  1 


3  2  U  (  x ,  y  ) 
3x2 


M-l  N-l 

y  y  c  (jTT/x  )2  sin  ( jirx/x  )  sin  ( kUy/y  ) 

>.  ik  max  max  m^x 

]  =  1  k= 1 


max 


(D.4) 


32U( x, y ) 
3y2 


M-l  N-l 

l  l  C  (ktr/y  )2  sin(jilx/x  )  sin(klTy/y  ) 
,  ,  ) k  max  max  max 

3*1  k= 1 


v  v 


By  substituting  equations  (D.4)  and  (D.5)  into  equation  (D.l)  the  following 
is  obtained 


O 


M-l 

N-l 

-  1 

l  [  (  j  1T/x 
k=l 

j  =  l 

Solving  for 

Clk 

r*  — 

jk  ' 

The  F.,  is 
]k 

computed  using 

(2.10),  and 

since 

M-l  N-l 

F(  x ,  y  ) 

=  l  l 

j  =1  k  =  l 

and  since. 

h  —  x  / ( M+ 1 ) 

x  max 

tr  - 

4  ( h  h  )  M 

_ X—X _ 

max  ]k 


max 


max 


-F(x.y)  ( D. 6 ) 


jk 


( D.  7  ) 


max 


jk~ 


max 


y  max 


max 


(D.2) 


.  .  .  J  /  F  (x,y)sm  rrijux/x  )sin(nkiry/y 

lk  (x  y  )  max  1  mav 

max  max  m=l  n=l 


max 


(  D.  8 ) 


By  using  equations  (D.8),  (D.7),  and  (D.3)  the  nodal  points  at  any  value 

of  U(x,y)  can  be  found  and  are  the  three  step  method  discussed  in 
Chapter  IV. 
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Appendix  E:  Mathematical  Explanation 
of  the  2-D  Complex  FFT 

To  show  that  only  the  imaginary  portion  of  the  FFT  array  must  be 
used  for  the  FFTCC  IMSL  routine,  the  following  mathematical  computation 
is  presented.  The  FFTCC  routine  computes  the  FFT  using  the  form 


A 

mk 


N-l 


l 


i2n(kn/N) 

A  e 
mn 


(E.  1) 


To  compute  the  two  dimensional  FFT  it  is  necessary  to  develop  two  trans¬ 
forms  similar  to  equation  (E.l)  using  different  notation  to  distinguish 
between  the  FFT  in  the  x-direction  and  the  FFT  in  the  y-direction.  These 
two  equations  are  A^1^  and  A*2'  and  defined  as  follows. 


A(*> 

mk 


N2-1 


l 

n  =  0 


i2rr(kn/N2) 

A  e 
mn 


(E.2) 


and 

A(2)  =  T  A  ei2*Um/Nl) 
jk  L  mk 

J  m=0 


(E.  3) 


where  N2  is  defined  as  the  number  of  transformed  points  in  the  x-direction 

(1) 


and  N 1  is  defined  as  the  number  of  transforms  in  the  y-direction.  A 


mk 


is  defined  as  the  Fourier  Transform  of  A  and  A^2^  is  defined  as  the 

mn  jk 


( 1) 


Fourier  Transform  of  A  and  the  double  Fourier  Transform  of  A  .  To 

mk  mn 


(2) 


solve  for  A  ,  one  can  substitute  equation  (E.2)  into  equation  (E.3)  and 
3  K 


obtain 


(E.4) 


© 


(2 

j 


Nl-1  /  N2- 1 

l  l 


m=0 


n  =  0 


A  e 
mn 


i2n  (  kn/N2  )  \  i2TT(jm/Nl) 


Using  Euler's  identity 


Nl-1 


l 


N2-1 

J  A  [  cos  (  2irkn/N2  )  +  isin  (  2irkn/N2  )  ] 
mn 


[cos(  2^3m/Nl )  +  isin  (  2TT3m/Nl )  ]  (E.5) 


By  multiplying  out  equation  (E.5),  the  following  is  obtained. 


Nl-1  N2- 1 

y  y  A  [cos(  2TTkn/N2  )cos  (  2it  jm/Nl ) 

„  „  nnn 

m=0  n=0 


+  isin(  2irkn/N2  )cos(  2tt  jm/Nl )  +  isin (  2tt  jm/Nl  )cos (  2irkn/N2  ) 

-  sin(  2irkn/N2  )sin(  2TTjm/Nl )  ]  (E.6) 

Note  that  the  only  term  required  to  compute  the  double  Fourier  sine 
transform  is  the  last  line  of  equation  (E.6),  -sin  (  2nkn/N2  )  sin  (  2tt  jm/N  1 )  . 
It  is  obvious  that  the  first  three  cosine  and  sine  terms  of  equation 
(E.6)  can  be  eliminated  by  deleting  the  real  portion  of  the  FFT  each 
time  it  is  computed  in  the  FORTRAN  program  FFT2D1  (See  Appendix  B). 
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